{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "sns.set_style(\"whitegrid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in log\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEYCAYAAADmugmLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8XGW5wPFftiZt03TfF1q6PKSFllKgQGVTdhBUQBZRwBUUr4joFRdE9N6LXr1XrrsiO7KIiFUWFRUQZGlDKS2dPqV0zdI2e9IkM8nMnPvHmZQ0ttmT95yZ5/v55NM5J2fmPJmeeZ953/MuWZ7nYYwxxgRRtusAjDHGmIOxJGWMMSawLEkZY4wJLEtSxhhjAsuSlDHGmMCyJGWMMSawLEkZY4wJLEtSxhhjAsuSlDHGmMDKdR1Ab7z++utefn7+vu1YLEbHbdM7nd+/5ubmqmXLlk10GNKA6Xyt9IaL6yps50zna8XKlf7r+B7291oJVZLKz8+nuLh433YkEtlv2/RO5/evpKRku8NwBlTna6U3XFxXYTtnOl8rVq70X8f3sL/XijX3GWOMCSxLUsYYYwLLkpQxxpjAsiRljDEmsCxJGWOMCaxQ9e4z6UdE7gTOA/ao6uEH+H0WcDtwDtAMXKWqrw1tlMYYV6wmZVy7Gziri9+fDcxP/XwS+OkQxGSMCQhLUhlqd0OUupaE6zBQ1eeBmi4OuQC4V1U9VX0ZGCMiU4cmuszU0pqgtLYZ3dVIyfYanttUybaqJtdhmRDY0xBlfVn9gL6mNfdlqOsfep18L8rxR7mOpFvTgZ0dtktT+yq6elIsFiMSifTphNFotM/P7auhPmci6bFxVyNPbXqFrbWtlDe0Ud0cp6o5wd7W5L8cf+i4Yfz4vTOGLD4TDo3RNl7ZUsMLm6v459tVbNq9l7ycLB784CEDdg5LUhmqIdpGYbbnOoxBYzNO/KtYPMFfI3v44xvlvLi5mvqWNgBGDsth7qRCFkwr4qTRBUwuKmBC4TAK8/MYmZ9DYX4ucycWMnbksH2vVVJSMqixmmBqjSdZs6OWFzdX8cLmKtaW1pNIehTkZXPM7HFceNQMTj1sEoma0gE7pyWpDJX0IDvLdRQ9UgbM7LA9I7XP9FBTLM59L2/nl89vobqplUmj8jlj4WTmjGjlvccvYvqY4WSH5GIwQyuZqnG3J6VXt9bQ0pYgOwsWzxjDtSfPZcW8CRx1yBjyc3P2PS/SVQN+L1mSylCe55EVjnJpJXCdiDwELAfqVbXLpj7zjuc2VXLTb9+gvD7KSQsm8tEVszlx/kRysrOIRCLMHDfCdYgmYHbWNO9LSi+9XU11UysAcyeO5INHz2DFvAksP3Q8o4fnDUk8lqQyVNLzyApAlhKRB4FTgAkiUgp8A8gDUNWfAU/idz/fjN8F/Wo3kYZLMunxvT8rP3n2beZNKuSRTx3PsXPGuQ7LBFBNUysvvV3NC5ureHFzFTtqmgGYNCqfkxdMZMW8CayYN4EpowucxGdJKkMlvWB07VTVy7r5vQd8ZojCSQuJpMfnH36dlWvLuezYmXzjvYsoyMvp/okmI7S0Jli1rWZfbWlDRQOeB4X5uRx36Hg+umI275o/gbkTCwPxRdaSVIZKhqe5z/SC53l87fF1rFxbzhfPFD59ytxAFDTGnXgiybqy+n1J6bXtdbQmkuTlZHHUrLHccNoCVsyfwOLpo8nNCcJX1/1ZkspQngdWdKWfu/+5jQdf3clnTp3LZ06d5zoc44DnebxduZcXN/tNeC+/XU1jLA7AwqlFXLViNivmTeCY2WMZMSz4KSD4EZpB4XecCN63JtN360rr+c8nI5xWPJkbzxDX4ZghtLshuq+m9OLmKnY3xACYOW445y2Zyop5Ezj+0PGMLwzfisOWpDJUiLqgmx6IJ5J88dG1TCjM53sXL7YmvjTXkBpE256YNu/ZC8DYEXmcMG8C75o3gRVzJzBrfPh7b1qSylB2Tyq93PvSdjbuauRnVyxjzIhh3T/BhE4snuDxNWU8tGonb3QYRHvsnPH7uoYXTylKuzFvTpOUiIwB7gAOBzzgo6r6ksuYMkUy6QWid5/pv4ZoGz94ZhMnLZjImYsmuw4nENKpbGmMtvHrV3bwqxe2sqcxxmFTRvHpU/xBtEtn7T+INh25rkndDjytqheJyDAg/HXTkIgnPXLS7BtXprr7xW00RON86UyxZr53hL5s2dMY5a4Xt3H/y9tpjMZZMW883//gEt41b0JG/T87S1IiMho4CbgKQFVbgVZX8WSahCWptLA3FudXL2zltOLJHD59tOtwAiHsZcvWqiZ+8fwWfvtaKfFEkrMPn8qnTj6UxTPGuA7NCZc1qTlAJXCXiCwBSoDPqepB1wToPLO1i9mq00WsLQ7JhL1/Iff4mjLqW9q49pS5rkMJkl6VLUEpV7QqyqPr63lxexO52VmcPq+QCxeNYVpRHjRWEImEZzawgXwPXSapXOAo4LOq+oqI3A58Gfj6wZ7QeWZrF7NVp42sHQzLy93v/bOZrcPF8zzuf3k7C6cWcdSszPyWfRC9Kltcliue5/H8W1X87Nm3eWlLNaMKcrn2lLlctWI2k0a5mYZoIHR8D/tbrrhMUqVAqaq+ktp+FP9CMkMgnkySY619ofbajjo27mrkP99/REbdo+iBUJQtz2zYzf/8ZRMbKhqYXJTPV845jMuOncWogqGZuDUsnHXwUtVdwE4RaR91+B5gg6t4Mk0i6ZFtBVuoPfZaKcPzcrjgyGmuQwmUoJctNU2tfPbBNXz83tVE4wm+e+Finv/SqXzypLmWoA7Ade++zwIPpHrfbMFmuB4yfscJ11GYvoonkjy9fhfvLp7EyHzXH+NACmTZ8sQbFdz8+/U0RNu44fQFXHvKXPLsg9glp1e3qr4OHO0yhkyUTHokPcixmlRovbylhuqmVt67eKrrUAIpaGVLZWOMm3+/nqfW72LxjNE8cNFyDptS5DqsULCvYBko4fnLxmfbF7jQemJdOSOH5XCKTHIdiumC53msXFvOLSvfpCmW4EtnCZ888dBAzjYeVJakMlAi6Scpq0mFk+d5/DWyh1MOm2TrRAVYazzJ1x9fz8Ord7J01hj++6LFzJs0ynVYoWNJKgPF4kkA8qx7XyhFKhrZ0xjjlAUTXYdiDqK+uY1rHyjhn29X89l3z+P60xbY4Pk+siSVgWLxBADDLEmF0rOb9gBwsiWpQNpR3czVd7/Kjppmvn/xEi5cNsN1SKFmSSoDxdpSNSn7ZhdKz2klC6cWMakovIM909W60nquvOtVkp7H/R9bzvJDx7sOKfTs7l0Gam/us5pU+DTF4pRsr+Ukq0UFzsZdDXz4zlcYMSyH3316hSWoAWI1qQzUavekQuu1HbXEkx7Hz7UCMEi2VO7lijtepSA3hwc/cRwzx4Vu0vXAsppUBmq/J2VJKnxWbaslOwubqy9A9jREueKOV/y5FD++3BLUALOaVAay5r7wWrW1huKpRTZ9TkDE4gmuub+E2uY2fnPN8cybVOg6pLRjNakM1GpJKpRa40nW7KzlmNnjXIdiUm5ZuYHXdtTxvYuX2Hpeg8SSVAaycVLhtL68nmhbkmPnWJIKgifeqODBV3dwzclzOdempxo0lqQy0L5xUtYFPVRe214LwNGHjHUcidndEOWrj69jyYzRfOGMBa7DSWuWpDJQNDVOaliuJakwWVdWz5SiAhsfFQDffiJCc2uC/7nkSJvFfJDZu5uB9kbbABiRZ//9YbKurN7uewTAq1tr+MPacq456VDmTrSOEoPNSqkM1NTqN/cNtyQVGo3RNrZWNXGEJSmnEkmPW1a+ybTRBVx7yjzX4WQEK6UyUGM0zrCcbOvdFyJvljfgebB4hiUpl377WikbKhq46Zxihg+zGeiHgiWpDLQ31kZhgQ2RC5P1ZfUA1tznUCLp8dNn32bRtCLOs958Q8aSVAZqiiUYmW/fAsNkXVk9U0cXMHFUvutQMtbT63extaqJz5w6jyxbi23IWJLKQI3ROIX5NmNBmKwvq2fRNKtFufTLf2zh0AkjOXPRFNehZBRr88lAe2NtFAaoJiUiZwG3AznAHap6W6ffzwLuAcakjvmyqj455IE6Eosn2FbdzDlHWBOTK7qrkdd31vG1c4tt8cIhZjWpDFTb1MbYEcNchwGAiOQAPwbOBhYCl4nIwk6HfQ14RFWXApcCPxnaKN3aUtlEIukxf7ItPe7Kw6t2kpeTxQeOsgUMh5olqQxU3dTK+MJgJCngWGCzqm5R1VbgIeCCTsd4QFHq8WigfAjjc27T7kYAxJKUE63xJL9bU8rpCyczbmRgPjcZw3lzX+qb9GqgTFXPcx1PuksmPWqbWxk/Mh9Iug4HYDqws8N2KbC80zG3AH8Wkc8CI4HTunvRWCxGJBLpU0DRaLTPz+2rrs75zzdryMmCWNVOIrUD19Tk4u8cSgNVtryytZra5jbed+T0gQvO9JjzJAV8DojwzjdlM4jqW9pIJL3UN8Ko63B66jLgblX9vogcD9wnIoer6kGzbH5+PsXFxX06WSQS6fNz+6qrc1a/upq5kwpZfHjnVtDBO2d3SkpKBjSWQTIgZcszG3ZTkJfNifNtNWQXnDb3icgM4FzgDpdxZJLqplaAIDX3lQEzO2zPSO3r6GPAIwCq+hJQAEwYkugCYNPuRrsf1UsDVbZ4nsdfNuzmXfMm2uBdR1zXpH4AfAno0SewcxNOujdXDIZ1u1oAaKreRXRcdhDev1XAfBGZg5+cLgUu73TMDuA9wN0iUoyfpCqHNEpHmlvj7Khp5qJldsO+l3pctnRVrmytiVFeH+XihYVB+KyExkCWzc6SlIicB+xR1RIROaUnz+nchOOiWSbstrRVABUsXTgf6sr2e/9cNOGoalxErgP+hN+9/E5VfVNEbgVWq+pK4AvAL0Xk8/idKK5SVW/Ig3Vg8569ACywmlSP9bZs6apceeXFrUAZF590BDPG2rLwPdXxPexvueKyJrUCOF9EzsH/ZlwkIver6hUOY0p7ZXXNAEwbM5zyOsfBpKTGPD3Zad/NHR5vwL9eMs5bu/0kNX+yzbbdCwNWtqzaVsv0McMtQTnkLEmp6k3ATQCpbzs3WoIafGW1LYwqyGX08LzM6scdUlurmsjJzmLWOCske2qgyhbP83hlaw0nzs+Y25+BZOOkMkxZXQvTxwx3HYbpoa3VTcwcO9wW1nOgtLaFqr0xltlKyE657jgBgKo+CzzrOIyMUFrbwoyxlqTCYmtlE7MnjHQdRmj1p2zZUNEAwKJpNjrGJft6lkE8z6Os1mpSYeF5Htuqm5hjScqJDeUNZGWBTLFOKy5ZksogexpjNMbiHGpLXofCnsYYza0JS1KORCoamDN+JCOGBaLBKWNZksog7XPAWU+xcNhS2QRgScqRyK4GiqdaU59rlqQyyL7uzJOs+SIMtlX7SWr2eEtSQy0WT1Ba28LcSfaFzjVLUhnkrT2NjB2Rx4TgTIlkurC1qolhudlMs3uIQ66stgXPw7r+B4AlqQyycZc/B5wtfR0OW6uaOGTcCFtkz4HtNf6g90PGW5JyzZJUhmiNJ3mzvIElM2wJ8rDYWmU9+1zZ2Z6krCblnCWpDBGpaKA1nmTpLBuYGAaJpMeO6mZLUo5sr26mIC+biaPyXYeS8SxJZYg1O2oBWDprjONITE/sbojSmkgyy5qbnNjVEGXa6OHWNB4AlqQyRMmOOiYX5TN1tN2ED4P25iab2NSNysYYE6wWFQiWpDJAMunx4uYqVsy1iTLDorTWX/drpk1h5URVY8ya+gLCklQG2FDRQE1TKycusCQVFjtrm8nKgumWpJyobIwxsdCSVBBYksoAz7/lL2K7Yp4lqbDYWdPC5FEF5OfakuVDLRpP0hiLW00qICxJZYC/RfZQPLWISaMKXIdiemhnbbPNVu9IXUsCwGpSAWFJKs2V17Wwenst5xw+xXUophfKaluYaWN0nGiIJQEYO9JmZgkCS1Jp7sl1FQCct2Sa40hMT7UlklTUt1inCUea2vwkVVRgs58HgSWpNPeHNypYNK3IBoWGSHldC0kPZlhNyommVr+5b1RBnuNIDFiSSmuRigbW7qzj/Uunuw7F9MI73c8tSbnQ1JqqSQ23mlQQWJJKYw+8sp1hudlctGyG61BML7wzkNea+1xoT1JWkwoGS1Jpam8szu9eK+O8xVMZM8JuAIfJztpmcrKzmDraemO60H5PqjDfalJBYEkqTT306g6aWhN8+LhDXIdiemlnTQvTxhSQm2MfTxeaWpOMys+1JVICwtlXBRGZCdwLTAY84BeqerureNJJtC3Bz5/fwglzx9us5yFUWtts96P6ob9lS3NrklHWsy8wXH5ViwNfUNWFwHHAZ0RkocN40sbDq3ZS2Rjjs++e7zoU0wc7a1vsflT/9KtsiSU8CobZTB9B4SxJqWqFqr6WetwIRADrhtZPTbE4P/77Zo6ZPZbjDh3nOhzTS7F4gsrGmM1+3g/9LVticY8Cm44qMAJRpxWR2cBS4JWujovFYkQikX3b0Wh0v20D966pYU9jjC+fOJ6NGzd2eay9f8Gzqz4KYJ0mBkhPypZ/KVfaEiSTCfts9MNAli3Ok5SIFAK/Ba5X1Yaujs3Pz6e4uHjfdiQS2W8705XVtfDYA9s4f8k0PnDS0m6P7/z+lZSUDGZ4pgfK6/wkNW2MNff1V0/Lls7lSuKpcsaOGmllSz90LFv6W6447T4kInn4F9EDqvqYy1jCzvM8vvH7NwH497MPcxyN6auKen8gr9Wk+qc/ZUss4ZGfZz0rg8LZ/4SIZAG/AiKq+j+u4kgXf3ijgmciu7nxDGG6fQsPrYp9zX32f9hX/S1b2hJJuycVIC6b+1YAHwbWicjrqX1fUdUnHcYUStV7Y9yy8k2WzBzDR981x3U4ph/K61oYOyKP4da7rD/6VbZYTSpYnCUpVX0BsNFy/ZRMetz4m7Xsjcb57oWLbQBiyFXUR60W1U/9LVtaE9a7L0js60LI3fHCFv6ulXz13GJkyijX4Zh+Kq/zZ5sw7rQmPAqsJhUYznv3mb4r2V7Ld59Wzlo0hY8c7376IxF5DP9ewFOqmuzF884CbgdygDtU9bYDHPNB4Bb8GQTWqurlAxJ0wFTURzlmto1vc6k14ZGfZzWpoLCvCyFVXtfCNfeXMHVMAd+5aDFZWYFo5vsJcDnwlojcJiLS3RNEJAf4MXA2sBC4rPPsACIyH7gJWKGqi4DrBzzyAGhujVPf0sZUq0k51ZbwGGbzJgaG1aRCqCkW5+P3rKalNcEDH1/O6OHBWFJAVZ8BnhGR0cBlqcc7gV8C96tq2wGediywWVW3AIjIQ8AFwIYOx3wC+LGq1qbOs2cQ/wxn9o2RsntSznieR9LD7u0GiCWpkEkkPa5/+HU27mrgzquOYcHkYN2HEpHxwBX4vavWAA8A7wKuBE45wFOmAzs7bJcCyzsdsyD12i/iNwneoqpPdxVH51kEesPFTBzRaJTX1ikArXW7iUS6HNc+YOe0WRX2l0h6AORakgoMS1Ihkkx63PTYG/xlw26+ef4iTpFJrkPaj4j8DhDgPuC9qlqR+tXDIrK6Hy+dC8zHT3IzgOdF5AhVrTvYEzrPItAbLmYyiUQi5BUVArs4fslhzByCpeP783em6+wkCc9PUtmWpALDklRIeJ7Ht57YwCOrS/m398znyhNmuw7pQH7ZeSyKiOSrakxVjz7Ic8qAmR22Z6T2dVQKvJJqLtwqIpvwk9aqAYo7EMrrW8jKgslFdk/KFatJBY/dHQyJ/33mLe56cRtXr5jN508L7BIc3z7Avpe6ec4qYL6IzBGRYcClwMpOxzxOqqlQRCbgN/9t6V+owVNRF2VCYT7Dcu1j6Up7krJ7UsFhNamA8zyP7/95Ez/6+2YuXjaDr5+7MCg9+fYRkSn495aGi8hS3hlIWQR02W6lqnERuQ74E/79pjtV9U0RuRVYraorU787Q0Q2AAngi6paPUh/jjPl9S1Mszn7nLIkFTyWpALM8zy+9ccId764lUuOnsl/fuCIoLaVnwlchd9U13GutEbgK909OdVE+GSnfTd3eOwBN6R+0lZFfZR5Ewtdh5HRLEkFjyWpgEokPb72+DoefHUnV50wm5vPWxjUBIWq3gPcIyIXqupvXccTRp7nUVHXwonzJ7gOJaNZkgoeS1IBFG1L8IXfrOWJNyr49Clz+eKZErgmvo5E5ApVvR+YLSL/UtuxWe6719SWpKk1YWOkHGvv3ZcT4M9bprEkFTC1Ta184t7VrN5ey01nH8anTp7rOqSeGJn619qq+qiyKQFgs004Fk9YTSpoLEkFyPbqJq66axVldS386PKlnLd4muuQekRVf57695uuYwmryqY4YOtIuWbNfcFjSSogXttRyyfuWU3C8/j1x5dzdIgmGRWR/+vq96r6b0MVS1hVpZKUzYDu1r7mPktSgdGjARkicnFP9pne8zyPX7+yg0t//jIj83N57NoTQpWgUkq6+THdqGyKk5OdxaRRmZWkgla2WE0qeHpak7oJ+E0P9pleiLYl+Nrj63m0pJSTFkzk9kuOZOzIYa7D6rVU7z7TD5XNcSaPys/EwjFQZYvNOBE8XSYpETkbOAeY3qlJpwiID2Zg6W5HdTPX3F/ChooG/u098/nce+aHtoASkR+o6vUi8gf89Z72o6rnOwgrVKqa4kwdkzn3o4JatrQnqWzr3RcY3dWkyoHVwPns32zTCHx+sIJKd3/fuIfrH34dz/O486qjefdhk12H1F/3pf79ntMoQqyyKc5RkzKqqS+QZcu+mlSOJamg6DJJqepaYK2IPKCqVnPqp0TS4/a/vsUP//YWh00p4udXLGPW+MGf7XqwqWpJ6t/nUvPvHYZfo1JVbXUaXAh4nkdVc4JpGVSTCmrZEreaVOB019z3iKp+EFgjIgdqxlk8aJGlmbrmVj730Os8t6mSC4+awbffdzjDh6XXEtUici7wM+Bt/Pn75ojIp1T1KbeRBVtNUyutCY+pGTRvX3DLFktSQdNdc9/nUv+eN9iBpLP1ZfVcc38Jexpi/Mf7D+fyY2cFegaJfvg+cKqqbgYQkbnAE4AlqS5U1Psr8mbYGCkrW0yPdNfcV5H6d/tgnFxEzgJux5/9+g5VvW0wzuPSI6t38rXH1zNh5DAeueZ4jpw5xnVIg6mxPUGlbMG/x2C6UF7XAmTWGCkrW0xP9agLuog08q+9turxb3x+QVV7vbaPiOQAPwZOx1/UbpWIrFTVDb19rSCKtiX45h/e5MFXd7Ji3nj+79KljC/Mdx3WoBCRD6QerhaRJ4FH8K+Xi0mzhQkHQ4bWpIDglS3evzQ8Gtd6Ok7qB/j/2b/Gv9dwKTAXeA24k9SCdL10LLC5/SIUkYeAC4DQJ6nS2mauvf811pXV8+lT5vKFMyS03ct76L0dHu8GTk49rgQyr+TtpYr6KLnZMD6EY+QGQCDLlvRsjQ+nniap81V1SYftX4jI66r67yLS7XpBBzEd2NlhuxRY3tUTYrEYkUhk33Y0Gt1vOwhKypr5zj/2kEh63HzqZI6f5bFJN7oO64AG6v1T1asHIJyMVVHfwoQRuYFdimWQOS9bOpYr2/b4tdodO3YQiVf18fRmIMvmniapZhH5IPBoavsiIJp6PGQV5Pz8fIqLi/dtRyKR/bZdSiY9fvLsZr7/110smDSKn314GXMmjOz+iQ51fv9KSvo3g5GIFAAfAxYB+26wqOpH+/XCaa6iLsqEkRk7jabzsqVjudI0vAYoZ9asWRTPnzgUp09LHcuW/pYrPZq7D/gQ8GFgT+rnw8AVIjIcuK6P5y4DZnbYnpHaFzr1LW188r7VfO/Pmzh/yTR+95kTAp+gBsl9wBT8lXqfw/8/tY4T3ahoaGFi5iapQJUtdksqeHr0yUi17b73IL9+oY/nXgXMF5E5+BfQpcDlfXwtZzbuauBT95VQVtvCN89fxEeOPyRdu5f3xDxVvVhELlDVe0Tk18A/XAcVZMmkx676KMdNK3IdihNBLVuyyNjPcOD0tHffDOCHwIrUrn8An1PV0r6eWFXjInId8Cf8bqJ3quqbfX09F/6wtpwvPfoGhQW5PPyp41h2SOhmLx9obal/60TkcGAXMMlhPIFX3dRKW8LL2JqUlS2mOz39ZNyF3/umfQr9K1L7Tu/PyVX1SeDJ/ryGC/FEku88vZFf/mMrRx8ylp986CgmFWXOGJcu/EJExgJfB1bir9T7dbchBVtFvT9GasKIzExSWNliutHTT8ZEVb2rw/bdInL9YAQUdNV7Y1z36zW8tKWaK48/hK+eu5BhuT29tZfeVPWO1MPngENdxhIW5XV+H4FMrUkRsLLFxkkFT08/GdUicgXwYGr7MqB6cEIKrjfL6/nEPaupbmrlexcv4aJlM1yHFCgiMh64Bb/pxsNvuvmWqmbctdJTu6wmFciyJXNvKwdPT6sAHwU+iH+PoQK/m+hVgxRTIP1lw24u/tlLeMCj15xgCerAHsLvoXUh/jVSBTzsNKKAq6iPMiw3m9EFGVsbz/iyxXStp737tuOv+7JPqkr+g8EIKkg8z+OOf2zlP5+KsHj6aH75kaPt/tPBTVXVb3XY/raIXOIsmhCoqI8ydXRBxvYIzeSyxfRMf9oYbiDNL6S2RJKbf7+eB1/dyTlHTOH7Fx+ZdstrDLA/i8il+HP3gf+t+E8O4wm8ivoWptiXns6clS2e3ZQKnP4kqbT+6tfcGuea+1/j+U2VXHfqPG44fUGmTlvTrQ6ThGYB1wP3p36VDewFbnQUWuCV10U5dk7GD13ozPkHzXkAZp/+JKm0/cpR19zK1XevYu3OOr5z4RFccsws1yEFmqqOch1DGCWTHrsbokwZXUAaf5z6wt4Ms093K/MeaBp98L9opOXs1hX1LXzkV6+yvaaZn16xjDMXTXEdUqiIyPnASanNZ1X1jy7jCbKqvTHiSY9powuAFtfhDKlMLFtM33S36GFGfUPeWdPMpb94mfqWNu65+liOnzvedUihIiK3AccAD6R2fU5EVqjqTQ7DCqz2daSmjB5OpiWpoJYtVoULnowdnNFZe4Jqao1g6VYUAAAZaklEQVTz0CeP4/Dpo12HFEbnAEeqahJARO4B1gCWpA6gfbaJqaMLoMFxMGZ/dlMqMDJ2cEZHpbV+gtobi3P/x5ZbguqfMR0e2xvZhfbZJqaOtt59xhxMxtekdjdEufQXL9MYbePXn7AaVD/9F7BGRP6O/130JODLbkMKrl0NUfJzsxk3chh7XAdjTEBldJJqiLZx5Z2vUtvUyoPWxNcvIpKFv7TCcfj3pQD+XVV39eC5ZwG3489YfYeq3naQ4y7EXxzvGFVdPSCBO1Re15LRA3mDyIZJBU/GJqlYPMGn7i1h85693HX1MSyeMab7J5mDUlVPRJ5U1SPwZ0DvERHJAX6MP+t1KbBKRFaq6oZOx40CPge8MoBhO7Wrvr37uQkaW08qODLynpTneXz5t+t4aUs1371oMSfaMtED5TUROab7w/ZzLLBZVbeoaiv+/H8XHOC4bwHf4Z2lxUOvoj7KtNHW29qYrmRkTepXL2zld2vK+MLpC/jAUTZR7ABajr/09zagCf++lKeqi7t4znRgZ4ft0tTr7CMiRwEzVfUJEfliTwKJxWJEIpFehP6OaDTa5+f2VCLpUVHfQl68iUgkMiTn7MzFOY3prYxLUv/cXMV/PbWRsxZN4bp3z3MdTro5c6BfUESygf+hlzNj5+fnU1xc3KdzRiKRPj+3p3Y3REl6Wzl87gyKiw8ZknN21p9zlpSUDHA0weDZSKnAyagktbshynUPrmHOhJF874NL7Ib1ABGRAuAaYB6wDviVqsZ7+PQyYGaH7Rmpfe1GAYcDz4oIwBRgpYicH+bOE+V1qTFSNrlsIFnREBwZk6SSSY8bf7OW5tY4j3zqeArzM+ZPHwr3AG34ixyeDSzE7+TQE6uA+SIyBz85XQpc3v5LVa0HJrRvi8izwI1hTlDwzmwTU8dYkjKmKxnTceKel7bxj7eq+Nq5C5k3qdB1OOlmoapeoao/x1+e48SePjFV47oOf0mPCPCIqr4pIrem5gFMS2W1fk1q+hjrOGFMVzKiOrGlci//9dRG3n3YJD603GY0HwRt7Q9UNZ5qlusxVX0SeLLTvpsPcuwpfYgvcMrqWijMz2X08DzXoZiO7JZU4DhJUiLy38B7gVbgbeBqVa0bjHN5nsfNv3+T/NxsbvvAEXYfanAsEZH22eeygOGp7fbefUXuQgum0toWpo8ZbtfjABuossX+V4LDVU3qL8BNqW/d38GfgPTfB+NEf3yjghc2V3HrBYts2fdBoqq2XHEvldW1MH2sNfUNgiErW8zQcJKkVPXPHTZfxr+PMeCaW+N8648bOHx6ER9afshgnMKYPimrbeboQ8a6DiPtDFXZYoZOEO5JfRR4uCcHdh6g2d1gxIfX1bKnMcaX3jWeTbqx34GmGxvM6UZDtI2GaNxqUoOvR2VLx3JlW4XfoWX79u0UxWza374ayLJl0JKUiDyDP6als6+q6u9Tx3wViPPOInld6jxAs6vBiPXNbTz28N94z2GTuOjkpb0NPyN0fv/SdYBm0FjPvv4Z6LKlY7lSk1cFVDB79myK54wbuKAzTMeypb/lyqAlKVU9ravfi8hVwHnAe1R1wPvU3PHCFhpjcW48s3c9zYwZbO1JaobVpPrEddlihpar3n1nAV8CTlbV5oF+/ebWOPe+tJ2zFk2heKp1LDPBUpaabcKa+wbeYJctZui5Gsz7I/zpbv4iIq+LyM8G8sUfLSmlvqWNj5946EC+rDEDoqyuhWG52UwYme86lHTUr7LF1pMKHle9+wZtZtdk0uPOF7aydNYYllnvKRNAZakxUtnZNhpnoA1U2WLD14Ij7aZFemlLNduqm7nqhNmuQzHmgErrWqzThDE9lHZJ6tGSUkYV5HLmogN1/jHGvfaalDGme2mVpPbG4jy9fhfnLZ5GQZ5NgmCCJ9qWoGpvzHr2BZStJxU8aZWk/vzmLlraElx41HTXoRhzQNazLxzsllRwpFWSeiaym8lF+Rw1yzpMmGAqtYG8xvRK2iSpWDzBc1rJe4onW68pE1g7avyhO4eMH+k4EmPCIW2S1MtbamhqTXBa8STXoRhzUDuqm8jPzWbSKBsjFUQ2Tip40iZJPaeV5Odmc8LcCd0fbIwj26ubmTVuhNX2A87GSQVH2iSpV7dVs3TWGOvVZwJtR00zh4wf4ToMY0IjLZJUQ7SNDeUNHDtnvOtQjDkoz/PYUdPMrHF2P8qYnkqLJFWyrZakB8fZ1PomwCr3xmhuTTBrnPXsCyq7JRU8aZGk1uyoJTsLjpw1xnUoxhzUjmrr2RcedlMqKNIiSW2oaODQiYWMGBaEhYaNObDtqSQ1y+5JGdNj6ZGkyhtYNM3WjTLBtr2mmawsW+wwyDzrgx44oU9StU2tlNdHWWiLG5qA21HdxLTRw8nPtR6oxvRU6JNUpKIBgIVWkzIBt73GHyNlgs/GSQVH6JPU21VNAMybVOg4EmO6tr3axkgZ01uhT1Lbq/xpZiaPKnAdijEHVdvUSk1TK3Mn2pepILM7UsET/iSVGsFv08yYINtStReAuZOs+3kYWGkSHOFPUtVNNu7EBN7be/xmaatJGdM7oU5SSc/z2/ntZrQJuLcr9zIsJ5sZY+1aNaY3Qp2kGqJJYvGkrXJqAu/tyr3MmTCSHGuWDja7KRU4TpOUiHxBRDwR6dP6GrUtcQAmWacJE3BbKpvsftQQ6m/ZkmV90APDWZISkZnAGcCOvr5GTUsCgElFtoCcCa7WeJLtNc12P2qIDETZYoLD5WR3/wt8Cfh9X19gX5KyVU5DTUTOAm4HcoA7VPW2Tr+/Afg4EAcqgY+q6vYhD7SPdtQ0kUh6HDrRalJDpN9liwkOJ0lKRC4AylR1rYj0+HmxWIxIJLJve09DFICa8m007wn17TUnotHofu+nCyKSA/wYOB0oBVaJyEpV3dDhsDXA0araLCLXAt8FLhn6aPtms/XsGzJ9KVs6lis7Sv1JgLdt3Ur+XruN0FcDWbYMWpISkWeAKQf41VeBr+BXx3slPz+f4uLifduNr1YxKj+XpUcs6nOcmSwSiez3fpaUlLgI41hgs6puARCRh4ALgH1JSlX/3uH4l4ErhjTCftq0u5GsLEtSA2Wgy5aO5UpF1m5gF3PmzKF4pi3901cdy5b+liuDlqRU9bQD7ReRI4A5QPs3nRnAayJyrKru6s05alsSTLSmvrCbDuzssF0KLO/i+I8BT3X3op1r3b0x0DXMVzbtZmphLju2vDVk5+yJINSk+2IoyhYTHEPe3Keq64BJ7dsisg2/Kaeqt6/VGEsyZoQlqUwhIlcARwMnd3ds51p3b3SuYfZX+RO7WTxrQpevOdDn7In+nNNRrbtLA1m2mOAI9Y2cptYkowryXIdh+qcMmNlhe0Zq335E5DT85pzzVTU2RLH1W3NrnG3VTRw2dZTrUEwP2HJSweN8KVtVnd3X5za1JSkabkkq5FYB80VkDn5yuhS4vOMBIrIU+DlwlqruGfoQ+27T7r14HhTbemdDrj9liw2TCo7Q16SKCpznWdMPqhoHrgP+BESAR1T1TRG5VUTOTx3230Ah8BsReV1EVjoKt9c2ptY7K55iScqYvghtCe95HntbE9bclwZU9UngyU77bu7w+IA3ysMgUtHAyGE5tmS8MX0U2ppULJ4knoSi4aHNsyYDRCoakSmjbCmZkLB7UsET2iTVEG0DoMhqUiagEkmP9eX1HDF9tOtQTC9l2YpSgRHaJNUY9SeXHWX3pExAbd6zl+bWBEtsUKgxfRbaJBVt8+ftK8jLcRyJMQe2dmcdgCUpY/ohxEkqCViSMsH1emkdowpymWMrR4eG3ZIKntAmqVjcr0nl54b2TzBpbu3OOpbMGGOdJkLIxkkFR2hL+Fjcr0lZkjJBFG1LsHFXI0tmWqcJY/ojtCV8rK09SVlznwme13fWkUh6HDVrrOtQjAm18Cap9ua+vND+CSaNvbylmqwsOHr2ONehmF7wbKBU4IS2hI9ZxwkTYC9vqWbRtCJG29ySxvRLeJOUdZwwARVtS/DajjqOmzPedSjGhF5oS3jrOGGCau3OOlrjSZYfaknKmP4KbQnfPpjXOk6YoHlxcxXZWXCs3Y8KHbsjFTyhTVKxeJLsLMjLsQENJlj+pns4atZYRo+w+1FhZeOkgiPUSSovJ4ssu5pMgOxuiLK+rIF3F0/q/mBjTLdCm6Ra40lyLUGZgPn7Rn/h4PccNtlxJMakh9AmqaTnkR3a6E26eiaym+ljhrNgcqHrUEwf2DCp4AltMZ9IetiUaCZI6ppbeW5TJecunmrN0CFn60kFR2iTVNLzyLaCwATIU+t30ZbwOH/JNNehGJM2nK0YKCKfBT4DJIAnVPVLvXm+1aRM0Dy+poy5E0eyaFqR61AyWn/LFhMsTmpSInIqcAGwRFUXAd/r7WskkpBjNSkTEFsq9/LK1hred+R0a+pzqP9li92UChpXzX3XArepagxAVff09gX85r4Bj8uYPrnv5e3k5WRxybEzXYeS6fpdtoCNkwoSV819C4ATReQ/gChwo6qu6u5JsViMSCQCQG1dHVl4+7ZN70WjUXv/BkBTLM6jq0s554ipTBpV4DqcTNfrsqVjuVJa2gTA1i1boC5/kENNXwNZtgxakhKRZ4ApB/jVV1PnHQccBxwDPCIih6pql3Xt/Px8iouLARi5poWc6tZ926b3IpHIfu9fSUmJw2jC656XttEYi3P1ijmuQ8kIA122dCxXticqgN3MOfRQiqfavcW+6li29LdcGbQkpaqnHex3InIt8FjqwnlVRJLABKCyp6+ftI4TJgAao2384vktnCoTOXLmGNfhZITBLFtsnFTwuLon9ThwKoCILACGAVW9eYFE0iPbspRx7OfPbaGuuY0bThfXoRhfv8sWsHtSQeLqntSdwJ0ish5oBa7srqmvs6QX4kFeJi28tbuRnz//Nh9YOp0jZox2HY7x9btsMcHiJEmpaitwRX9ew3r3GZfaEkm+9Ns3GJmfy1fPtfuiQTEQZYsJFmeDefvLH8xrWcq4cdtTG1mzo44fXraU8YXWCyxdWJUreELbYmYTzBpX7n5xK796YStXHn8I77UpkNKSzd0XHCGvSbmOwmQSz/O468Vt3PrHDZyxcDJfO2+h65CMSXshT1KWpczQaIrF+eYf3uSR1aWcvnAyP7x8KXk5VpU3ZrCFNkklPQ9bOd4MtkTS4+n1u/iPJzZQXh/lulPnccPpC2z4Q5qycVLBE9okZTWp9CEiZwG3AznAHap6W6ff5wP3AsuAauASVd02mDHtbojyxBsV3P/ydrZUNSGTR/Hby5ey7JBxg3laExBWtARHeJOUB7l2IYWeiOQAPwZOB0qBVSKyUlU3dDjsY0Ctqs4TkUuB7wCXDFQMDdE2tlc1s7mykZLttbyou9hWtwXPg8UzRvOjy5dy9uFTybHakzFDLrRJyqZFShvHAptVdQuAiDyEv9RCxyR1AXBL6vGjwI9EJKungzTL6lp46NUd1DW3Ud/yzk9DSxs1za3UNbftO3bksBwWjB/GDact4OwjpjJvki0Db4xLoU1SC6cWkRVrdB2G6b/pwM4O26XA8oMdo6pxEakHxtPFdDcdZ7b+69uN/PSflYzIy6ZwWDaFw3IozM9mRmE2xeOHM6VwFFNH5TG9KI+Zo/Noa41RUBCnrXonkeqB/FMPzsWM9DYL/r+aM2Ekk0bmMmmUjX0LitAmqe9ctNg+YOagOs5sXVwMnznX6/FihJ1nhx8KYTtnus6Yv3BaEfdcNIsxI4a5DsWkWB9a41oZ0HGlwBmpfQc8RkRygdH4HSh6zFbLNSacQluTMmljFTBfRObgJ6NLgcs7HbMSuBJ4CbgI+JtNGmpMZrCalHFKVePAdcCfgAjwiKq+KSK3isj5qcN+BYwXkc3ADcCX3URrjBlqVpMyzqnqk8CTnfbd3OFxFLh4qOMyxrhnNSljjDGBZUnKGGNMYFmSMsYYE1iWpIwxxgSWJSljjDGBleWFaG76kpKSSmC76zjS2CHLli2b6DqIgWDXyqCza8X0VL+ulVAlKWOMMZnFmvuMMcYEliUpY4wxgWVJyhhjTGBZkjLGGBNYlqSMMcYEliUpY4wxgRXKWdBF5CzgdiAHuENVb3McUqiIyJ3AecAeVT3cdTxBJCK3AJ8AKlO7vpKarX0wzjXk17OIbAMagQQQV9WjB/ucYWBly4EdqMwQkXHAw8BsYBvwQVWtFZEs/PfwHKAZuEpVX0s950rga6mX/baq3tPduUNXkxKRHODHwNnAQuAyEVnoNqrQuRs4y3UQIfC/qnpk6mewEpTL6/nU1N9mCQrn/xdBdzf/WmZ8Gfirqs4H/so767ydDcxP/XwS+CnsS2rfAJYDxwLfEJGx3Z04dEkK/4/brKpbVLUVeAi4wHFMoaKqzwM1ruMwgF3PQWL/FwdxkDLjAqC9JnQP8L4O++9VVU9VXwbGiMhU4EzgL6pao6q1wF/owZflMCap6cDODtulqX3GDLTrROQNEbmzJ9/4+sjV9ewBfxaREhH55BCcLwysbOmdyapakXq8C5icenyw97FP728o70kZMxBE5BlgygF+9VX8Jopv4Rfm3wK+D3x06KIbdO9S1TIRmQT8RUQ2pr4tG9NrquqJyKDMsRfGJFUGzOywPSO1z5heUdXTenKciPwS+OMgheHkelbVstS/e0Tkd/hNXZmepKxs6Z3dIjJVVStSzXl7UvsP9j6WAad02v9sdycJY3PfKmC+iMwRkWHApcBKxzGZNJP60LV7P7B+kE415NeziIwUkVHtj4EzGLy/L0ysbOmdlcCVqcdXAr/vsP8jIpIlIscB9almwT8BZ4jI2FTz+RmpfV0KXZJS1ThwHf4fFwEeUdU33UYVLiLyIPCS/1BKReRjrmMKoO+KyDoReQM4Ffj8YJzE0fU8GXhBRNYCrwJPqOrTg3zOwLOy5eAOUmbcBpwuIm8Bp6W2AZ4EtgCbgV8CnwZQ1Rr8pvNVqZ9bU/u6ZEt1GGOMCazQ1aSMMcZkDktSxhhjAsuSlDHGmMCyJGWMMSawLEkZY4wJrIxMUiKSEJHXRWS9iPxGREb08vl7e3n83SJy0QH2Hy0i/5d6fJWI/Cj1+BoR+UiH/dN6cz4TDL29Tg7yGtNE5NHU4yNF5Jz+R2aCRkRmiMjvReQtEXlbRG5PjdXqfNyzIvIvEwJ3LD/STUYmKaAlNfvz4UArcE3HX6YGoQ36e6Oqq1X13w6w/2eqem9q8yrAklSGUtVyVW3/gnMk/vIHJo2klrZ4DHg8NaP4AqAQ+A+ngQVEGKdFGmj/ABaLyGz8QXyvAMuAc0TkBOArQBb+gMd/b3+SiPwv/ojpXcClqlopIp/An5p+GP5Atg+ranPqKaeJyJeBIuAGVf2jiJwC3Kiq53UMKLWW0V78NVqOBh4QkRb8OeU+oarvSx13OvBpVX3/wL4lZrCkrrM7gQn4a1Vdrao7RGQu8AAwEn/k/vWqWpg6/o/AUcCtwHAReRfwX6r6sIM/wQy8dwNRVb0LQFUTIvJ5YGuqLPgVsATYCAxvf5KIXA3cBNQBa4FYav/F+EtiJPBnezhp6P6UgZepNSkARCQXf+2Tdald84GfqOoioA34Dv4FdCRwjIi0T0U/ElidOu45/AsC4DFVPUZVl+CPWO84k8Ns/PnRzgV+JiIF3cWnqo8Cq4EPqeqR+CO5DxORialDrsYv8Ex4/BC4R1UX4yel/0vtvx24XVWPwJ8dej+ppSNuBh5OtQJYgkofi4CSjjtUtQHYAXwBaFbVYvxyZhnsm7brm8AK4F3461+1uxk4M1UOnT/o0Q+yTE1Sw0XkdfwEsAP/mwrA9tT6JwDHAM+qamVqupQHgPZvJEn8FSkB7se/SAAOF5F/iMg64EP4F1+7R1Q1qapv4U8Zclhvg1ZVD7gPuEJExgDHA0/19nWMU8cDv049vo93rp3jgd+kHv+685NMxjoFv4xBVd8A3kjtX8475VMr75RHAC8Cd6dadnKGMNZBkanNfS2pmsk+IgLQ1MfXa59b6m7gfaq6VkSuYv8ZfzvPP9XX+ajuAv4ARIHfpBKoMSa8NgD7dawSkSJgFn6TcK+o6jUishy/1aZERJapavWAROpAptakeuJV4GQRmZBaVvoy/KY98N+39ovqcuCF1ONRQIWI5OHXpDq6WESyU/ceDgW0h3E0pl4X8G+kA+XA1/ATlgmXf+LPrg3+NfKP1OOXgQtTjy/t/KSU/a4Fkzb+Cozo0KM3B3/9sruBp/HLGETkcGBx6jmv4JdP41PlzcXtLyYic1X1FVW9GT/JdVw2I3QsSR1Eamr5LwN/x78pWaKq7VPRNwHHish6/HtWt6b2fx3/4nkR/yZnRzvwE99TwDWqGu1hKHfj38N6XUTab5o+AOxU1Uiv/zAzlEakZoxu/7kB+CxwdWp29Q8Dn0sdez1wQ2r/PKD+AK/3d2Bh6lq4ZCj+ADP4Us3478f/IvsWsAm/peQr+ItvFopIBL+cKUk9pwK4BX9m8hfx74G3++/UDP7r8b8UrR2iP2VQ2CzoIZQaD7FGVX/V7cEmFFJj9VpSK5xeClymqhe4jssY1zL1nlRoiUgJfk3uC65jMQNqGfCj1JiZOtJrqXpj+sxqUsYYYwLL7kkZY4wJLEtSxhhjAsuSlDHGmMCyJGWMMSawLEkZY4wJrP8HwNlb3kZSO60AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prob = np.arange(0, 1.0, 0.001)\n",
    "odds = prob / (1-prob)\n",
    "logits = np.log(odds)\n",
    "\n",
    "f, (ax1, ax2, ax3) = plt.subplots(1, 3)\n",
    "ax1.plot(prob, logits)\n",
    "ax1.set_xlabel('Probability')\n",
    "ax1.set_ylabel('Logit')\n",
    "#ax1.set_title('L vs. P(Y)')\n",
    "\n",
    "ax2.plot(logits, prob)\n",
    "ax2.set_xlabel('Logit')\n",
    "ax2.set_ylabel('Probability')\n",
    "#ax2.set_title('P(Y) vs. L')\n",
    "\n",
    "ax3.plot(odds, logits)\n",
    "ax3.set_xlabel('Odds')\n",
    "ax3.set_ylabel('Logit')\n",
    "#ax3.set_title('L vs. odds')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logits[int(len(logits)/2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prob=0.9500000000000001 logit=2.944438979166442\n",
      "prob=0.995 logit=5.293304824724491\n",
      "prob=0.999 logit=6.906754778648553\n"
     ]
    }
   ],
   "source": [
    "x=len(logits)-50\n",
    "print('prob={} logit={}'.format(prob[x], logits[x]))\n",
    "x=len(logits)-5\n",
    "print('prob={} logit={}'.format(prob[x], logits[x]))\n",
    "x=len(logits)-1\n",
    "print('prob={} logit={}'.format(prob[x], logits[x]))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
